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SUMMARY 

A noniterative finite-difference numerical method is presented for the solution of the incompressible 
Navier-Stokes equations with second-order accuracy in time and space. Explicit treatment of convection 
and diffusion terms and implicit treatment of the pressure gradient give a single pressure Poisson equa- 
tion when the discretized momentum and continuity equations are combined. A pressure boundary condi- 
tion is not needed on solid boundaries in the staggered mesh system. The solution of the pressure Poisson 
equation is obtained directly by Gaussian elimination. This method is tested on flow problems in a 
driven cavity and a curved duct. 


INTRODUCTION 

A finite-difference numerical method is presented in this paper for time-dependent solution of the 
incompressible Navier-Stokes equations. Owing to the special role of the pressure which results in the 
absence of a time derivative term of the pressure in the continuity equation, calculation procedures for 
incompressible flow differ significantly from that for compressible flow. The continuity equation in the 
form of a constraint prohibits straightforward integration of the system of flow equations with time. To 
overcome this difficulty Harlow and Welch (ref. 1) have presented the Poisson equation method for pres- 
sure by taking a divergence of the momentum equation and combining it with continuity. This has been 
extended to three space dimensions by Williams (ref. 2). 

Chorin (ref. 3) has introduced a fractional step method, which is now more often called the pro- 
jection method. This scheme has brought forth a variant known as the velocity-pressure method (refs. 4 
and 5) and it has become very popular for fluid flow and heat transfer problems in two or three space 
dimensions. This method generates an intermediate velocity field directly from the momentum equation 
alone, which does not satisfy the continuity constraint. The intermediate velocity is corrected succes- 
sively until the new velocity field becomes asymptotically divergence-free. The solution procedure is 
iterative in nature. Kim and Moin (ref. 6) introduced a pressure-related Poisson equation as a combina- 
tion of references 1 and 3. 

A direct coupling of the momentum and continuity was also proposed by Moin and Kim (ref. 7). 
This solves the momentum and continuity equations simultaneously to obtain the flow variables, say 
velocity components and pressure, at the advanced time level. The artificial compressibility method 
(ref. 8), developed for steady flows, has been modified to yield time-dependent solutions (refs. 9 and 10). 

In this paper we present a direct coupling method on a staggered grid, which also calls for the solu- 
tion of the discrete pressure Poisson equation. Time integration uses an Adams-Bashforth explicit 
approximation for convection and diffusion terms, and a Crank- Nicolson implicit approximation for the 
pressure term. The present method is different from the fractional step method since it uses no inter- 
mediate variable to be corrected at a later stage. 



FORMULATION 


The system of equations for incompressible viscous fluid flow can be written in dimensionless form 


as 


^ + V • (uu) = — Vp + — V 2 u (1) 

dt Re 

V • u = 0 (2) 


where u is the velocity, p the static pressure, t the time, and Re is the Reynolds number. Equa- 
tion (1) is the momentum equation in which each velocity component has its own time derivative. In 
other words, the momentum equation is written in a time evolution form. The continuity equation (2), 
which is a constraint as opposed to an evolution equation, should be satisfied at any time. 

First, we start with a finite-difference approximation of the time derivatives only. Let us express 
the finite-difference approximation of the system (eqs. (1) and (2)) as 


u n +l __ u n 

At 


C(u) — D(u) G(p) — 0 


V • u 


n+l _ 


(3) 

(4) 


where C(u), D(u), and G(p) are the time approximations of the convection, diffusion, and pressure 
gradient terms, respectively, and the superscript n is the time level. There are various ways to describe 
these terms in order to maintain second order accuracy in time. If the Crank-Nicolson approximation is 
used for all terms 


C(u) = l|v • (uu) n+1 + V • (uu) n ] 

(5) 

D(u) = _L(v 2 u n+1 + V 2 u n ) 
Re 

(6) 

G(p) = !(v P n+1 + Vp n ) 
2 

(7) 


This is a nonlinear simultaneous system of equations, hence, iteration using a block inversion or sequen- 
tial scalar inversion is required. To avoid the nonlinearity the convective term V • (uu ) n+1 in C(u) can 
be approximated by 


V • (uu) 


n+l 


dx. 
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where Au r = u^ +1 — u£. Cartesian tensor notation is invoked in the above expression, with subscripts 
p, q, and r taking values 1,2,3, and with the summation convention applying. The system (eqs. (3) 
and (4)) with the use of equation (8) can be written in component form as 


where 
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This system requires block matrix inversion at each grid point. Furthermore, the coefficients of the 
matrix should be calculated at each time step because of the implicit treatment of the convection terms. 
Consequently, the direct solution procedure must be repeated every time level. If convection is treated 
explicitly using the Adams-Bashforth approximation, equation (5) becomes 


C(u) 


1 

2 L 


3V • (uu) n - V 


(uu)”- 1 ] 


(10) 


Then, the operators M and N in equation (9) become 


1 - At(2Re) _1 V 


In this case, the matrix 


coefficients are to be determined once and remain unchanged throughout the entire calculation. However, 
a simultaneous solution for u, v, and p is needed, which requires a block matrix inversion. 


In the approach which will be used in this paper, the diffusion as well as the convection terms are 
treated explicitly using the Adams-Bashforth method. The diffusion term becomes 



Then, equation (9) becomes 
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here, and are a-components of C(u) and D(u) described in equations (10) and (11). The 
operators M and N in this case become the identity operator, so that a simultaneous solution proce- 
dure for u, v, and p can be avoided if the momentum equation is combined with the continuity equa- 
tion in an appropriate manner. The purpose of this paper is to utilize the time approximation in 
equation (12), with convection and diffusion terms explicitly treated and the pressure implicitly 
approximated. 


For convenience, the space discretization will be described in a two-dimensional Cartesian 
coordinate system, but it is readily extendable to three-dimensional space. Figure 1 shows the staggered 
mesh system for the spatial finite-difference formulation. In the staggered grid, the discretization of 
equation (12) yields 


n+l n _n+l _n+l 
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(15) 


where 


(Ru)y = 


i+lj 
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- C x (ij) + D x (ij) 


( Rv )ij 




2Ay 


- C y (ij) + D y (ij) 


Details of the central difference formulation for the convection and diffusion terms, which are C (ij) and 
D a (ij), can be found in reference 11. In order to combine the discrete momentum equation with the con- 
tinuity equation (15), we rewrite equations (13) and (14) at mesh points (i— 1 j) and (ij—1), respectively, 
to give 
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into 


n+1 ~ n+1 ^ n+1 
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n n+1 „ n+1 ^ nH 
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where C w = C e = Ax“ 2 , C 3 = C n = Ay~ 2 , C c = -2(Ax~ 2 + Ay“ 2 ), and 

l [( Ru )ij - ( Ru )i-ij + ( Rv )ij - ( Rv )ij-ll 


(16) 


( R P)ij = ( 2At ) 


The entire system of equations (eqs. (13) to (15)) now reduces to a single scalar equation (16), which is 
the discrete Poisson equation for pressure. New flow variables at (n+1) At can be obtained algebraically 
by putting the pressure solution of equation (16) into equations (13) and (14). 


The discrete Poisson equation (eq. (16)) in a five-point stencil form is solved by direct inversion 
using Gaussian elimination. A pivoting strategy is not needed because the matrix is diagonally dominant. 
The storage required in this inversion is approximately 2 X min(II,JJ) X min(II,JJ) X max(II,JJ), where 
II, JJ are the number of cells in the x and y directions, respectively. In this formulation the pressure 
boundary condition is not needed on the boundary such as a wall where the velocity is described. This 
should not really be a surprise, since pressure is introduced into the incompressible formulation only for 
the purpose of linking the momentum and continuity equations in their primitive variable form. The 
pressure boundary condition issue arises when it is needed to solve the differential pressure Poisson equa- 
tion, which is written to be V 2 p = d 2 u q v r /dx q dx r . In the present approach where the momentum and 
continuity equations are directly combined in their discrete form, the boundary condition is already taken 
care of by the momentum equation. However, the pressure boundary condition can be obtained by a 
comparison of the directly coupled system and the differential Poisson equation. The boundary condition 
is equivalent to dp/dn. — 0 for an inviscid case and dp/dn = Re” 1 ^^ /Sn 2 for a viscous case, where 
n is the are normal variable to the wall. However, the matrix equation (eq. (16)) is singular since all of 
the discrete equations are not linearly independent. This can be elucidated by looking at continuity in a 
discrete integral form over cell surfaces. The following expression can be derived by a mass conservation 
principle, 


II JJ 

E E 

i+1 j = l 
i j *m 


(u • ASJij = -(u ■ AS)^ 


(17) 


where AS is the area vector with outward normal to the cell surface. The above expression says that any 
one of the discrete continuity equations can be constructed by adding up all of the continuity difference 
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equations over the rest of the cells. Therefore, one of the equations should be discarded and the pressure 
value at the relevent cell can be set arbitrarily. This is a generic feature of incompressible flow calcu- 
lations, in which the pressure itself is to be adjusted around the arbitrarily fixed value to render the 
appropriate pressure gradient. The above equation (eq. (17)) is exact in our central difference formulation 
on the staggered grid. Also, equation (17) holds exactly for nonuniform grids and in generalized coordi- 
nates if the discrete continuity is written using any conservative scheme. F or global mass conservation, 
in which equation (17) holds exactly, the integral form of the continuity equation is recommended for a 
stretched grid or for coordinate systems other than Cartesian. This can be done easily by discretizing a 
surface integral of the mass over an individual cell. In the present work we set C w = C g = C g = C n 
= Rp = 0 and C c — 1 at (II, JJ) to fix the value p(II,JJ) to be zero. 

Because of the explicit treatment of the convection and diffusion terms the method presented here 
is subject to a stringent numerical stability limitation. This requires use of a time step within the 
stability limit. We choose the time step to meet the convection and diffusion limits, which have been 
referred to as Courant and Neumann limits, respectively. The Courant number C r and Neumann dif- 
fusion number D n are defined as 



r \ 


D = At 
n Re 
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c r = At 

| u | ^ | v | 

^Ax Ax^ 

and 

1 + 
,Ax 2 

> _ 

<< h “‘ 

to 


In the present algorithm using three time levels, the upper limit of D n is 1/4. However, the value of C r 
will be chosen smaller than unity in a heuristic way. 


EXAMPLE PROBLEMS 

The first example is for unsteady flow inside a square cavity after the upper wall is set in motion 
impulsively. As a second example, fluid flow in a curved duct of a square cross section is presented. The 
flow is assumed fully developed, which eliminates downstream variable dependency. 


Driven Cavity Flow 

Driven cavity flow has been preferrably chosen to test the numerical scheme for incompressible 
Navier-Stokes equations mainly because of the following two reasons. Firstly, since velocity is clearly 
described on the flow boundary, no ad hoc boundary condition on the velocity is necessary, in contrast to 
channel flows in which in-and-out flow boundary conditions can be specified in many different ways. 
Secondly, there is no primary flow direction. Any a priori assumption regarding directivity of the flow, 
which often allows one-dimensional approximation, is not appropriate to apply in this recirculating flow. 
Therefore, the finite-difference scheme should be more general than the simple upwind differencing 
method based on one-dimensional analyses. The majority of driven cavity flows have been emphasised for 
the steady case. Unsteady computations can be found in literatures (refs. 10, 12, and 13) for flows driven 
by implusively started or oscillatory motion of the lid. 

The present example demonstrates flow development with time when the upper lid, y = 1, is set in 
a uniform motion impulsively. Two cases are studied with the Reynolds numbers of 5000 and 10 000, 
which are reasonably large for incompressible laminar flow calculations. Strict central difference formu- 
lation is maintained in both diffusion and convection. No numerical device such as an artificial damping 
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term or switching of the difference formulation of the convection terms is applied in this study. A 120 by 
120 uniform grid is used throughout. For high Reynolds number flow the Courant number limit is more 
stringent than the Neumann limit. Therefore, only convection stability is considered taking C r = 0,6, in 
which u and v are set to be unity and zero, respectively. This yields At = 0.005. This value is not 
an upper limit and the time step could be increased on the same grid without affecting numerical 
stability. 

Bruneau and Jouron (ref. 14) report from their steady calculation that there would be no steady 
solution for a high Reynolds number flow. They explain that the solution does not converge and if con- 
verges, the solution becomes dependent on the grid size for Reynolds numbers larger than 5000. The cal- 
culation shows that the appearance and then vanishing of minor vortices repeats in the corner region of 
the lower part of the cavity during iterations. Goodrich et al. (ref. 15) show that the high Reynolds 
number flow inside the cavity is unsteady in its time asymptotic behavior. To achive this asymtotic 
oscillatory flow a long term calculation should be performed after an impulsive start. Several thousand 
nondimensional time units are required in their computation. The objective of the present study is to 
apply the current numerical method to a physical problem, rather than infer the ultimate oscillatory flow 
phenomena. Therefore, we focus our attention only on the early development of the impulsively started 
driven-cavity flow for nondimensional time units of 210 for Re = 5000 and 300 for Re = 10 000. 

The velocity, length, and time are made dimensionless by the lid velocity U , side dimension L, 
and L/U Q , respectively. The Reynolds number Re is defined to be U Q L/V, where v is kinematic vis- 
cosity. Immediately after impulsive movement of the lid at uniform velocity u = 1, an enormously steep 
velocity gradient forms in the vicinity of the lid, which gives rise to a huge drag force. The drag coef- 
ficient C d is defined as 


C, 


1 f 1 

r \ 

du 

Re jo 



dx 


(18) 


The drag coefficient appears to decrease precipitously to a steady value in a short period of time (this is 
not shown here but a similar behavior can be found in ref. 10 for Re — 400), which easily misleads one to 
conclude that the flow has reached a steady state. However, it takes a large amount of time to reach the 
steady state, especially for high Reynolds number flow. The steady state would be stationary, in which 
no change occurs with time, or oscillatory. Figure 2 shows the velocity profiles of u and v at the 
geometric centers, x = 0.5 and y = 0.5, at the end of present calculation. At t = 210 for Re = 5000 
the present results agree very well with the steady solution of Ghia et al. (ref. 16). The local maxima of 
u and v near the bottom and right walls are —0.418 at y = 0.079 and —0.536 at x = 0.954; whereas 
theirs are —0.436 at y = 0.07 and —0.554 at x = 0.953, respectively. The comparison for Re = 10 000 
at t = 300 in figure 3 shows a quantitative difference. The local maxima of u and v at the geometric 
centers of the cavity are —0.407 at y = 0.062 and —0.501 at x = 0.963; whereas theirs are —0.427 at 
y = 0.055 and —0.543 at x = 0.969, respectively. This indicates that further time integration is needed 
to reach the steady solution. As will be shown later, the local flow still changes at t = 300, especially in 
the corner regions. 

Figure 4 shows streamfunction contours at six different times for Re = 5000. The vortex center is 
located just underneath the lid immediately after the implusive motion and moves towards the right wall 
as shown at t = 1. Then, the vortex bounces against the wall and heads downward. The contour at 
t = 5 shows the downward movement of the vortex with its outer contour distorted. The distortion of 
the outer streamlines allows a bulge of attached vortex at the right wall, which is not shown here but it 
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occurs between t = 5 and 10. The vortex, which forms near the middle of the right wall, moves down 
and resides at the lower right corner as the secondary vortex. As time proceeds the secondary vortices 
form progressively in the lower-left and upper- left corners. The streamline contours at t = 100 and 210 
show little difference. 

The increments of the streamfuntion contour A tjj are 0.0015 at t = 1, 0.004 at 5, and 0.005 at 
later times. The contours in the corner regions where the secondary vortices reside are made with a dif- 
ferent scale since the magnitude of contour level of the secondary vortex is much lower than that of 
primary vortex. A Aifr of 2x10 4 is used at t — 15 and 5xl0~ 4 afterward. Moffatt eddies (ref. 17) 
are observed in the lower right corner for times greater than 30. The Moffatt eddies constitute a vortex 
cascade which consists of a sequence of eddies of decreasing size and rapidly decreasing intensity occurring 
near a corner between plane boundaries at rest. The innermost vortex at the corner will hereafter be 
called the tertiary vortex. The tertiary vortex is very small and weak for Re = 5000, so no attempt is 
made to put more contours in it using a finer scale. 

At Re = 10 000 the flow undergoes much more complicated processes. Figure 5 illustrates the 
time evolution of the flow using streamfunction contours. The early flow development, t — 1,3, and 5, is 
the same as for Re = 5000. At t = 5 the vortex is seen bulging from the right wall, which later forms 
the secondary vortex at the lower right corner after t = 20. The detatched vortex at t = 10 appears 
relatively strong and it rolls on the bottom wall (not shown), the left wall as shown at t = 20, and then 
finally stays at the upper left corner. This rolling vortex becomes the secondary vortex in the upper 
corner. At t = 30 the contour plot shows that the primary, secondary, and tertiary vortices are well 
identified at three corners of the cavity. However, the flow is far from steady and even the primary vor- 
tex changes its shape and the center long after t = 30. This can be seen in all of the contour plots in 
figure 5. As time procedes further the primary vortex appears to be stabilized qualitatively as shown at 
t = 100 and 300. 

As in the case of Re = 5000, different scales are used in the streamfunction contour plots of the 
primary and secondary vortex and no further smaller scale contour is added in the area occupied by ter- 
tiary vortex. The tertiary vortices appear in two corners of the lower part, which are much larger in size 
than that for Re = 5000. As shown at times larger than 90 the secondary and tertiary vortices near the 
lower right corner seem to become stable as time elapses. However, the vortices in the lower left and the 
upper corner do not tend to be steady. They change incessantly. They may take on oscillatory behavior 
in the long run. 

After a series of spiral motions the primary vortex center finds its ultimate position at (0.5166, 
0.5333) for both Re = 5000 and 10 000. The present calculation indicates that the primary vortex 
center does not change in either case after about t = 120. This agrees excellently with the result of Ghia 
et al. which shows the vortex centers are (0.5117, 0.5352) and (0.5117, 0.5333) at Re = 5000 and 10 000, 
respectively. The locus of the primary vortex center is shown in figure 6 for Re = 10 000. (No attempt 
is made to interpolate adjacent values of streamfunction to yield the exact vortex center. Instead, the 
maximum streamfuntion value is searched only on the grid point positions.) The vortex center is traced 
from t = 0.5 with a time increments of 0.5. 


Fluid Flow in a Curved Duct 


Fluid flow in a curved duct has been studied in many engineering applications both for developing 
and for fully-developed cases. Most of the applications are for the steady flow. Unsteady cases have been 
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also reported for a pulsatile fluid flow in a curved pipe (refs. 18 and 19), which is relevant to the blood 
flow through the ascending aorta. 


In the present work the flow is at rest initially and it is assumed that a sudden axial pressure 
gradient is imposed. Also a fully-developed flow assumption is made. This kind of fluid flow could be 
observed far downstream of a curved duct, whose entrance is connected to a reservior, and an exit valve 
is opened suddenly. Because of the fully-developed condition postulated the derivatives with respect to 
the downstream variable p are to be neglected except the pressure. The cylindrical coordinates (x, y, 
and <p ), along with the corresponding velocity components (u, v, and w), are illustrated in figure 7. The 
curvature ratio 6 is defined to be a/R, where a is the gap, r Q - r i? and R is the radius of curvature. 
The nondimensional Navier-Stokes equations are given in cylindrical coordinates as 


where B = 1 + 6x and 
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(19) 

(20) 
( 21 ) 

( 22 ) 


The Rq is the Reynolds number based on a reference velocity W r because for fully-developed case 

the mean velocity W m is only obtained as a part of the solution. The reference velocity has been chosen 

as 


W, 


f y/2 

aG 


where 


G = - 1 dp' 
R d<p 


where p is the density, p' the dimensional static pressure, G is the axial pressure gradient which is the 

driving force maintaining the flow in the duct, and the pressure p is nondimensionalized as P /(pW r j. 

The above system of equations is relevent to a rotating fluid flow with an axisymmetric geometry and 
heat transfer problems because the rotating velocity component and the energy are formally decoupled 
from the pressure gradient. 
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In the fully-developed flow case only the cross sectional components of the momentum equation 
(i.e., u and v equations) are to be combined with the continuity to yield a discrete Poisson equation for 
the pressure. Due to the area change in the radial direction the continuity is derived to be the surface 
integral of the mass over the cell rather than using equation (22) directly. This, as mentioned before, 
exactly satisfies the global mass conservation in a discrete sense. The discrete Poisson equation is con- 
structed as before, using the Adams-Bashforth approximation for the convection and diffusion and 
Crank-Nicolson for the pressure, as described in section 2. The axial momentum equation (eq. (21)) is 
integrated using the Peaceman and Rachford type ADI as 
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(24) 


where the quantities with the superscripts c and a denote the values at (n+l/2)At by the Crank- 
Nicolson and the Adams-Bashforth approximations, respectively. Since u and v at (n+l)At are 
already obtained by equations (19) and (20) and the surface integration version of equation (22), it is 
trivial to obtain these at (n+l/2)At by the Crank-Nicolson method. 

The Reynolds numbers of 120 and 200 are considered with a 60 by 40 grid in the cross plane (x,y). 
The curvature ratio 6 is taken to be 1/6.45. The time step At is more restrictive due to the 
decoupling effect between the w-momentum equation and the pressure (hence the continuity), and is 
chosen to be 0.0002. Numerical solution is sought over the entire domain of the cross-section, 

—0.5 < x,y < 0.5. Computations for both R Q = 120 and 200 are performed up to t = 24, which is con- 
sidered sufficient for a steady state. 

Calculations reveal that the flows for both Reynolds numbers are symmetric about the plane 
y = 0, which will be called the plane of symmetry hereafter. Figure 8 illustrates the main and the 
secondary flow development for R Q = 120. The contour plots are made using only the half domain of the 
cross-section. The upper half is for the streamfunction and the bottom half is for the w velocity. The 
cross-flow streamfuntion ^ is defined as 
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It is well known that the velocity maximum for steady flow occurs near the outer bend due to a 
centrifugal effect. Surprisingly, quite the reverse is observed in an early stage of flow evolution, such that 
the w maximum occurs near the inner bend as shown at t = 1 and 2 in Figure 8. This maximum 
migrates towards the outer wall as the flow develops further, which is illustrated at t > 4. It is observed 
that the w maximum, which is initially located on the plane of symmetry (shown at t = 1,2,4 in the 
plot), shifts from it and forms double peaks in the entire cross-section. This is shown for t > 8. The 
axial flow boundary layer is readily recognized along the outer and bottom (and top) walls. The w dis- 
tribution near the inner wall is not of boundary layer type because higher velocity fluid particles are 
squeezed into the outer bend due to the centrifugal force. The contour plots are made with different 
scales as time changes since the flow variable values are very small in the beginning. The values of Aw 
are 0.1 at t — 1, 0.2 at t — 2, and 0.25 at later times. 

Immediately after the sudden imposition of the pressure gradient, the secondary flow already takes 
on a vortex-pair pattern in such a way that fluid particles travel toward the inner bend along the bottom 
and top walls of the duct, collide and separate at the middle of the inner bend, and reenter through the 
plane of symmetry toward the outer wall. The strength of the secondary flow is negligible in the very 
beginning. The streamlines in the top half of the cross-section in figure 8 show the secondary flow 
history. (A^ is taken to be 10 4 at t — 1, 10~ 3 at 2, and 0.005 at t > 4.) The secondary flow for 
R o = 120 maintains a single vortex-pair from the beginning with the intensity increasing. It is worth 
noticing that a minor vortex-pair appears near the outer bend at t = 14 and 20. This additional small 
vortex-pair, which is of Taylor-Gortler type, remains minor at this R Q . 

For the case of R Q = 200 the flow exhibits more complexity. Figure 9 shows the axial velocity 
contours over the entire cross-section with Aw = 0.1 at t = 2 and 0.25 at the other times. As in the 
previous case, the w maximum occurs in the vicinity of the inner bend and gradually shifts toward the 
outer bend. Migration of the w maximum forms a boundary layer near the outer wall as shown at 
t = 6 and the high velocity contour becomes distorted severely. This can be seen in figure 9 for t > 8 as 
an eruption of the wall boundary layer from the middle of the outer wall. This axial flow complexity is 
closely associated with the secondary flow behavior, which will be seen later in the secondary flow con- 
tour plot. The w velocity profiles at t — 24 are drawn in figure 10. For R Q = 120 the migration of 
the maximum velocity toward the outer wall can be seen on the plane of symmetry, while the w max- 
imum for R q — 200 is pushed toward the center. Again, the velocity maximum for R Q = 200 is seen to 
shift toward the outer bend and off the plane of symmetry, (y is about 0.14 in fig. 10.) 

The secondary flow is complex for R Q — 200. The contour plots in the upper half of the 
cross-section in figure 11 are for the streamfunction. The increment A^ is taken to be 0.001 at t = 2 
and 0.005 at t > 4. One vortex pair, which is maintained in earlier times (shown at t = 2 and 4), 
becomes two vortex-pair. The additional vortex near the outer wall is not minor, but as strong as the 
primary vortex. The same scale of A if: is used for both the primary and the additional vortices. The 
strong vortex near the middle of the outer bend is considered responsible for the appearance of the w 
maximum near x = 0 on the plane of symmetry since it conveys fluid particles away from the outer wall 
through the symmetry plane. 

The v velocity contour is shown in the bottom half of the cross-section. Four extrema (for the 
full cross-section) in the v velocity correspond to the one-vortex-pair structure in the secondary flow as 
shown at t — 2 and 4. Likewise, eight extrema of the v are seen clearly corresponding to the two- 



vortex pair for t > 8 in figure 11. The Av is 0.006 at t = 2 and 0.03 at t > 4. The contour of v at 
R Q = 120, which is not shown here, exhibits four large valleys and peaks and additional extrema due to 
minor vortices at the outer bend are insignificant. 

A relation between the reference velocity W r and the mean velocity W m can be found as 

W m = QW r where Q = w dx dy 

The flowrate Q is based on W r . The Reynolds number defined as Re = W m a/V becomes Re = QR q . 
The Re is more conveniently used in experiments because the W m is a direct measure. The same is 
true in the analysis of spatially developing flow when an inlet velocity profile, hence the Re, is given. 
Figure 12(a) is for the flowrate Q which is also interpreted as the time history of Re. The values of Q 
at t = 24 are 2.8012 and 3.4645 for R Q = 120 and 200, respectively, and the Reynolds numbers become 
Re = 336.14 and 692.9 accordingly. These values are confirmed by the steady solutions by the artificial 
compressibility method mentioned earlier. Figure 12(b) presents the pressure difference between the 

outer and the inner wall on the plane of symmetry, which is |p Q — p. /(,W r 2 ). The dip for R Q = 200 

appears as a response to the change in the secondary flow structure from one-vortex pair to two-vortex 
pair. The flowrate history for R 0 = 200 in figure 12(a) also bears the vestige of the same response. 


CONCLUSIONS 

With the explicit treatment of the convection and the diffusion and the implicit treatment of the 
pressure, the direct coupling of the discrete momentum and the continuity equations renders a single pres- 
sure equation. Direct inversion of the matrix by the Gaussian elimination provides efficiency and speed of 
the computation since iteration is not needed each time step. The present method is tested on the driven 
cavity flow of high Reynolds numbers and the curved duct flow to give second order accurate solutions in 
time and space. Details of the secondary and tertiary vortices are captured in the driven cavity. The 
solution for Re = 5000 agrees very well with the published data. For Re = 10 000 the result does not 
show the steady behavior at least within the time units used, but qualitativly agrees with the steady 
data. It is found that the velocity maximum is located near the inner bend in the early stage of flow 
development. Time evolution of the main and secondary flow towards the steady state and changes in 
vortex pattern are clearly simulated. 
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Figure 5. — Streamf unction contours versus time at Re = 1 0 000. 
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Figure 12. — Time history of flowrate and pressure difference in a curved duct. 
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